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1  Executive  Summary 

This  final  report  summarizes  the  research  carried  out  from  May  2010  to  October  2012.  The 
tasks  outlined  in  the  original  proposal  for  two  years  are:  (1)  literature  review  and  development  of 
a  peridynamic  solver,  (2)  development  of  an  automatic  volume  calculation  routine  using  the 
partition  of  unity  principle,  (3)  numerical  testing  for  different  grid  widths  to  horizon  ratios,  (4) 
development  of  an  approach  to  add  another  material  variable  in  the  given  approach,  (5) 
preparation  of  the  imaging  system  for  experimental  testing,  (6)  development  of  peridynamic  and 
finite  element  coupling  module,  (7)  parallelization  of  the  solvers,  (8)  verification  of  the  coupled 
code,  and  (9)  preparation  of  the  imaging  system.  As  for  the  literature  review  for  the  peridynamic 
formulation,  intensive  study  has  been  conducted.  A  peridynamic  solver  has  been  developed  in 
Fortran,  and  the  effect  of  the  different  grid  sizes  and  horizons  are  investigated  for  brittle 
materials.  The  numerical  methodology  has  been  extended  to  accommodate  the  plasticity 
behavior  of  materials.  The  peridynamic  code  is  implemented  using  a  graphics  processing  unit 
(GPU)  for  highly  parallelized  computation.  A  comparison  of  stresses  and  strains  by  finite 
element  analysis  (FEA)  and  peridynamic  solutions  is  performed  for  a  ductile  material.  A 
multiscale  procedure  is  proposed  to  bridge  the  material  models  by  means  of  peridynamic  bonds 
and  the  corresponding  macroscale  material  model  used  in  the  finite  element  analysis.  A 
numerical  scheme  to  couple  the  peridynamic  nodes  and  finite  elements  is  developed  in  Matlab 
and  is  verified.  From  the  point  of  view  of  experimental  work,  a  high  speed  camera  setup  was 
prepared,  and  impact  tests  were  conducted.  The  digital  image  correlation  technique  (DIC)  was 
selected  to  analyze  the  high  speed  response  of  the  material.  During  the  phases  I  and  II,  four 
research  papers  have  been  submitted  to  international  journals.  Three  papers  have  been  published, 
and  the  last  one  is  under  review. 

2  Introduction 

Fracture  is  one  of  major  concerns  in  the  engineering  field  for  a  long  time,  and  people  have  made 
substantial  efforts  in  order  to  understand  material  failure.  Inglis  and  Griffith  made  contributions 
to  the  early  development  of  fracture  analyses,  and  Irwin  extended  the  Griffith  approach  by 
developing  the  energy  release  rate  [1],  Fracture  mechanics  provides  a  more  reliable  methodology 
for  engineering  designs  than  the  traditional  strength  based  approach.  However,  classical  fracture 


mechanics  has  its  limitations.  For  example,  a  pre-existing  crack  needs  to  be  defined,  and  the 
fracture  propagation  zone  is  required  to  be  small  compared  to  geometrical  dimensions  [2,  3], 

In  order  to  simulate  crack  growth  and  material  damage,  various  numerical  methods  have  been 
studied.  The  cohesive  crack  model  addresses  relationships  between  the  crack  opening 
displacement  and  the  cohesive  tractions  resisting  the  separation  of  cracks  [4-6],  The  cohesive 
zone  models  have  been  incorporated  into  finite  element  (FE)  models  using  interface  elements 
and  contact  surfaces.  The  partition  of  unity  finite  element  method  (PUFEM)  was  presented  by 
Melenk  and  Babuska  [7].  Belytschko  and  collaborators  [8,  9]  extensively  investigated  the 
partition  of  unity  principle  for  simulations  of  fracture  problems  and  developed  the  extended 
finite  element  method  (XFEM).  The  XFEM  allows  the  discontinuity  not  constrained  to  element 
boundaries,  and  it  can  model  the  discontinuity  without  remeshing  [10].  Mariani  and  Perego  [11] 
presented  a  method  for  the  simulation  of  quasi-static  cohesive  crack  propagation  using  the 
XFEM  for  quasi-brittle  materials.  Cox  [10]  proposed  enrichment  functions  to  represent  the 
discontinuity  using  an  analytical  investigation  of  the  cohesive  crack  problem.  Considering  the 
minimization  of  the  total  energy  of  the  global  system,  Meschke  and  Dumstorff  [12]  proposed  a 
variational  form  of  XFEM  for  the  propagation  of  cohesive  and  cohesionless  cracks  in 
quasi-brittle  solids.  Meshffee  and  particle  methods  [13]  demonstrate  capability  for  numerical 
simulations  of  material  failures.  Molecular  dynamics,  for  example,  is  capable  of  investigating 
nonlinearities  in  the  vicinity  of  cracks,  the  bond  breaks  between  atoms,  and  the  formation  of 
extended  defects  [14-16].  Research  on  coupling  finite  element  methods  and  meshless  methods 
[17,  18]  has  also  been  conducted. 

A  great  deal  of  research  effort  has  been  made  to  study  many  fracture  problems.  One  common 
benchmark  problem  characterized  by  the  mixed  mode  fracture  is  the  test  of  a 
double-edge-notched  concrete  specimen  conducted  by  Nooru-Mohamed  et  al.  [19],  The  test  of 
Nooru-Mohamed  was  adopted  by  De  Borst  [20]  in  the  discussion  of  computational  modeling  of 
concrete  fracture.  For  the  analyses,  the  finite  element  smeared-crack  approach  with  the  gradient 
Rankine  plasticity  model,  Cruch-Crack  model,  and  Ottosen’s  model  has  been  used  for  numerical 
studies  of  Nooru-Mohamed’ s  experiment  by  Di  Prisco  et  al.  [21],  A  comparative  study  of 
three-dimensional  constitutive  models  for  the  double-edge-notched  test  was  performed  by 


Pivonka  et  al.  [22],  Gasser  and  Holzapfel  [23]  employed  the  cohesive  crack  model  with  the 
PUFEM  for  the  numerical  modeling  of  the  test.  The  XFEM  was  utilized  by  Cox  [10],  Meschke 
and  Dumstorff  [12],  and  Unger  et  al.  [24]  for  the  simulations.  An  adaptive  mesh  refinement 
technique  applied  to  a  nonlocal  version  of  anisotropic  damage  model  was  employed  by  Patzak 
and  Jirasek  [25],  Rethore  et  al.  [26]  used  a  hybrid  analytical  and  XFEM  to  study  the  propagation 
of  curved  cracks  in  the  double-edge-notched  concrete  specimen. 

As  a  reformulated  theory  of  continuum  mechanics,  peridynamics  eliminates  the  spatial 
derivatives,  and  it  is  valid  regardless  of  discontinuities  [27,28],  Therefore,  peridynamics  is  useful 
to  solve  problems  involving  spontaneously  emerged  discontinuities.  With  the  general 
applicability  of  peridynamics,  many  applications  of  the  method  have  been  studied.  Silling  and 
Askari  [29]  wrote  the  first  paper  on  the  numerical  simulation  using  the  peridynamic  model.  In 
their  work,  the  bond-based  peridynamics  is  employed  to  study  the  convergence  in  a  fracture 
problem  and  impact  of  a  sphere  on  a  brittle  target.  Dayal  and  Bhattacharya  [30]  studied  the 
kinetics  of  phase  transformations  using  peridynamics  without  any  additional  kinetic  relation  or 
the  nucleation  criterion.  By  adding  pairwise  peridynamic  moments,  Gerstle  et  al.  [31]  proposed  a 
micropolar  peridynamic  model.  Other  applications  of  peridynamics  include  modeling  of  the 
structural  responses  under  extreme  loading  [32],  structuralstability  and  failure  analyses  [33], 

fracture  analyses  [34-37],  and  nucleation  analyses  of  a  crack  in  a  solid  body  [38],  Peridynamics 

also  has  been  applied  to  study  the  effect  of  fiber  directions  in  composites  on  the  fracture  and 
damage  evolution  considering  anisotropic  properties[39-43],  A  generalized  formulation  of 
bond-based  peridynamics  was  introduced  by  Silling  et  al.  [44],  which  is  called  the  state-based 
peridynamics.  The  convergence  of  peridynamic  states  to  classical  elasticity  was  studied  by 
Silling  and  Lehoucq  [45],  and  it  is  shown  that  the  peridynamic  stress  tensor  converges  to  the 
Piola-Kirchhoff  stress  tensor  as  the  length  scale  goes  to  zero.  Using  the  state-based  peridynamic 
method,  Warren  et  al.  [46]  studied  the  elastic  deformation  and  fracture  of  a  bar.  Littlewood  [47] 
presented  fragmentation  of  an  expanding  tube  modeled  with  state-based  peridynamics.  Foster 
et  al.  [48,  49]  studied  viscoplasticity  and  the  failure  criterion  for  peridynamic  states. 


Compared  with  FEM,  peridynamics  is  computationally  expensive.  Macek  and  Silling  [50] 


implemented  peridynamics  in  a  commercial  finite  element  analysis  code,  ABAQUS,  using  truss 
elements.  The  conventional  FE  mesh  is  coupled  with  the  peridynamic  truss  mesh  using  the 
embedded  element  feature  available  in  the  finite  element  analysis  code.  Lall  et  al.  [51]  used  the 
peridynamics  based  finite  element  model  to  study  shock  and  vibration  reliability  of  electronics. 
Kilic  and  Madenci  [52]  presented  a  coupling  approach  using  overlapping  regions  in  which  both 
peridynamic  and  FE  equations  are  utilized.  Agwai  et  al.  [53]  and  Oterkus  [54]  employed  the 
submodeling  approach  to  couple  the  FEM  with  peridynamics.  In  their  approach,  the  global 
analysis  by  means  of  finite  element  analysis  is  performed  first,  and  then  peridynamics  is  used  for 
submodeling.  A  morphing  strategy  based  on  the  energy  equivalence  was  proposed  by  Lubineau 
et  al.  [55], 

In  this  research,  we  introduce  a  coupling  approach  of  discretized  peridynamics  with  finite 
elements.  Different  from  the  approach  in  [50,  51]  implementing  peridynamic  model  in  the 
framework  of  the  conventional  FEM  and  the  submodeling  approach  [53,  54],  the  peridynamic 
subregion  is  directly  coupled  to  the  finite  element  subregion  in  the  present  approach.  An 
interface  element  is  introduced  to  calculate  coupling  forces  instead  of  using  overlapping  regions 
[52]  or  the  morphing  strategy  [55]  to  couple  peridynamic  and  FE  subregions.  Depending  on  how 
coupling  forces  are  divided  to  FE  nodes  of  an  interface  element,  we  further  discuss  two  types  of 
coupling  schemes. 


Figure  1 :  Schematic  of  peridynamics. 


3  Theory 

3. 1  Peridynamic  theory 

In  the  peridynamic  theory,  the  equation  of  motion  of  a  material  point  at  x  in  the  reference 
configuration  at  time  t ,  as  shown  in  Figure  1,  is  written  as  [27,  56,  57] 

/?u(x,0  =  f  f(Tj,^)dVx,  +b(x,0,  (1) 

JHx 

where  p  is  the  mass  density,  u  is  the  displacement  vector,  f  is  a  pairwise  force  vector  that 
the  material  point  at  x'  exerts  on  the  material  point  at  x,  Hx  is  a  neighborhood  of  the 


material  point  at  x ,  and  b  is  the  body  force  density  field.  The  relative  position  vector  in  the 
reference  configuration  is  expressed  as  [29] 

<r=x'-x,  (2) 

and  the  relative  displacement  vector  at  time  t  is  written  as 

77  =  u(x',0-u(x,0.  (3) 


For  each  material,  a  scalar  5 ,  called  the  horizon,  is  assumed  to  exist  to  determine  the  interacting 
spatial  range  between  the  material  point  at  x  and  the  material  point  at  x'  such  that 

f(77,^)  =  OV^  ifP£B >S,  (4) 

where  P-  P  is  the  Euclidean  norm.  The  pairwise  force  vector  f  in  the  bond-based 
peridynamics  is  expressed  as 


f  0?,#)  =  /(*,£); 


V  +  4 


Ptj  +  g  P’ 

where  /  is  a  scalar-valued  pairwise  force.  The  bond  stretch  s  is  defined  as 

P/7  +  £P-P£P 


s(t,TJ,^)  =  - 


P£P 


(5) 

(6) 


where  P£P  is  the  original  bond  length  in  the  reference  configuration,  and  P^  +  ^P  is  the 
deformed  bond  length.  If  the  bond  stretch  5  =  0,  then  there  is  no  pairwise  force  /  between 
material  points.  Figure  2  shows  a  microbrittle  material  model  defined  for  peridynamic  bonds. 
The  bond  force,  which  is  a  scalar  function  of  the  bond  stretch  s ,  is 

/(*7»  £)  =  cs(h  £)/*(?,  h,  4), 


(7) 


where  c  is  the  micromodulus,  and  //  is 
the  scalar  function  to  determine  the  bond 
failure.  The  micromodulus  c  can  be 
determined  by  equating  the  strain  energy 
density  in  the  classical  elasticity  [29]  or  by 
numerical  calculations  [56].  The 
numerically  determined  micromoduli, 
obtained  by  considering  a  finite  number  of 
bonds  within  the  horizon,  are  summarized 


Figure  2:  Microbrittle  material  model. 


in  Tables  1  and  2  for  one-  and  three-dimensional  models,  respectively.  The  critical  stretch  for 
bond  failure  is  denoted  as  s0 .  Once  a  bond  fails,  it  cannot  sustain  force  any  more.  The  critical 

bond  stretch  s0  for  microbrittle  materials  is  obtained  by  setting  the  work  required  to  break  all 
the  bonds  per  unit  fracture  area  identical  to  the  energy  release  rate  Gf  [29]: 


*0=' 


'9  kS' 


(8) 


By  considering  broken  bonds,  damage  dependencies  can  be  introduced  into  the  critical  bond 
stretch  [29,  34],  The  scalar  function  /u  is  related  with  the  critical  bond  stretch  as 


\ 1  if  s(t',  tj,  £ )  <  s0  for  all  0  <  t'  <  t, 
1 0  otherwise. 


(9) 


To  solve  the  peridynamic  equation  of  motion,  the  material  domain  is  discretized  with  a  number 
of  nodes.  The  distances  of  two  adjacent  nodes  are  identical  over  the  domain  and  denoted  as  Ax  . 
Therefore,  the  volume  representation  of  each  node  is  (Ax)3.  The  peridynamic  equation  of 
motion  after  discretization  is  written  as 


N  u 


M  =  Zf(V’£)F/+b/’ 


(10) 


j= i 


where  iij  is  the  acceleration  of  the  node  /  at  time  t ,  is  the  pairwise  force,  is 

the  total  number  of  nodes  within  the  horizon  of  the  node  I ,  and  b)  is  the  body  force  at  time  t . 
The  damage  index  of  a  node  is  defined  as  [29] 


(11) 


(p(xI,t)  =  l-J^r - 

H/ 

Ay* 

j= i 

In  order  to  consider  the  volume  reduction  of  a  node  that  has  an  intersection  with  the  horizon 
boundary,  a  volume  reduction  scheme  [58]  is  introduced  as  follows: 


Horizon  8 

Micromodulus  C}  (N/m  6 ) 

2Ax 

E 

2Ax4 

3Ax 

2  E 

9Ax4 

4Ax 

E 

8Ajc4 

5Ax 

2  E 

25Ax4 

Table  1 :  Micromodulus  c, 

in  one-dimensional  domain. 

Horizon  8 

Micromodulus  C3  (N/m 6 ) 

2Ax 

0.302942  Ea 

Ax4 

3Ax 

0.052385  Ea 

Ax4 

4Ax 

0.017290  Ea 

Ax4 

5Ax 

0.006819  Ea 

Ax4 

Table  2:  Micromodulus  c3  in  three-dimensional  domain. 
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(12) 


where  (S  -  r )  is  the  distance  from  which  the  volume  is  reduced,  and  r  is  set  to  be  the  half  of 
the  grid  spacing  Ax  in  the  numerical  implementation. 


3.2  Finite  element  formulations 

The  equation  of  motion  in  the  conventional  continuum  mechanics  is  derived  from  the  principle 
of  linear  momentum.  The  temporal  change  rate  of  linear  momentum  is  equal  to  the  force  applied 
on  the  body  as  [59] 

PUt  =  Fjj  +  fB  inQ>  (13) 

where  p  is  the  density,  ut  is  the  acceleration  field,  er  is  the  Cauchy  stress  tensor,  and  fB 
is  the  body  force  density  field.  The  essential  boundary  condition  rw  and  natural  boundary 
condition  Ty  are  defined,  respectively,  as  [60] 


ii 

on 

(14) 

W  =  Fi' 

on  r/5 

(15) 

where  the  surface  of  the  body  r  =  r„ur/>  r«nr/  =  0 ,  and  Hj  means  the  components  of  the 
unit  outer  normal  vector  on  F  .  The  constitutive  equation  for  continuum  is  stated  as 

~  ^ijk-Fkn  (16) 

where  Cijkl  is  the  elastic  constitutive  coefficient,  and  the  components  of  strains  are  defined  as 

£ij=\{uij+uj.i)-  (17) 

Applying  the  principle  of  virtual  work  to  Equation  (13),  we  have 

=  (18) 


where  Sll  is  the  virtual  displacement.  After  integrating  by  parts  and  applying  the  divergence 


theorem,  the  weak  formulation  is  obtained  as  [61] 

-f  (JnSuj  ,<fQ+  f  <j..n/SujdS  +  f  f^Su.dO.  -  f  piidudQ.  =  0.  (19) 

Jn  y  Jr  v  1  '  Jn  '  ’  Jn  11 

Considering  the  symmetry  of  the  stress  tensor  (cr.  =  a  ji)  and  applying  the  boundary  conditions, 
we  have 

f  pii.8u.dO.  +  [  a,8s,dQ  =  f  Fis8uidS+  f  fB8u,d Q.  (20) 

Jn  '  '  Jn  u  v  h f  1  1  Jn  '  1 

Substituting  the  constitutive  law  in  Equation  (16)  into  Equation  (20),  we  obtain  the  finite  element 
formulation 

\dM,Su,dn+\cIJklsuSsIJdn  =  \  F’su,ds+y,‘su,dn.  (21) 

In  the  finite  element  analysis,  the  displacements  within  each  element  are  interpolated  by  means 
of  shape  functions  as  [60] 

u(e)  =H(e)U(e),  (22) 

where  H(e)  is  the  displacement  interpolation  matrix,  the  superscript  e  represents  the  element 
e ,  and  the  nodal  displacement  vector  U(e)  is  expressed  as  U,c-)r  =  {w,  r,  w,  •  •  •  un  vn  wn\  for  an 
element  of  n  nodes.  The  strain  vector  is  evaluated  by 

eie)  =B(e)U(e),  (23) 

where  B(e)  is  the  strain-displacement  matrix  which  is  written  as 
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Substituting  Equations  (22)  and  (23)  into  Equation  (21),  we  have  the  weak  formulation  in  matrix 
form  as 


MU  +  KU  =  F,  (25) 

where  M  is  the  mass  matrix,  K  is  the  stiffness  matrix,  and  F  is  the  force  vector.  The 


assembled  matrices  following  the  convention  of  direct  stiffness  method  [60]  are  summarized  as 


M  =  ^M(e),  M<e)  =  J  (e)/?(e)H(e)rH(e)JQ(e),  (26) 

e 

K  =  ^K(e),  K(e)  =  Jn(e)B(e)rC(e)B(e)JQ(e),  (27) 

e 

F  =  2X(e)+ZFae)’ 

e .  .  (28) 

Fs(e)  =  Jr(e)H(e)rF5(e: 'dS(e) ,  Fge)  =  Jn(e)H(e)rf  ^ B(e)dn(e) . 


4  Coupling  between  the  peridynamic  and  finite  element  subregions 


4. 1  Coupling  schemes 

To  gain  the  efficiency  from  finite  element  analyses  and  exploit  the  generality  of  peridynamics  in 
the  presence  of  discontinuities,  a  domain  is  partitioned  into  a  conventional  FE  subregion  and  a 
peridynamic  subregion  as  shown  in  Figure  3.  With  the  coupling  approach,  the  subregion  where 
failure  is  expected  can  be  modeled  using  peridynamics.  The  overall  computational  cost  is 


reduced  by  using  FEM  for  the  remainder  of  the  domain. 
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Figure  3:  Partitions  of  the  domain.  The  FE  subregion  and  the  peridynamic  (PD)  subregion  are  bridged  by 
interface  elements. 


In  the  conventional  FE  subregion,  the  lumped  mass  matrix  is  formed  by  distributing  the  total 


element  mass  to  nodes  of  the  element  [60],  The  internal  forces  on  FE  nodes  can  be  calculated  as 

p/n>  =  yyint(e)  =  ^K(e)U(e),  (29) 

e  e 

where  K(e)  is  the  element  stiffness  matrix.  The  equation  of  motion  of  an  FE  node  is  obtained  as 

M/U/=F/art-F“,  (30) 

where  M;  is  the  lumped  mass  of  node  I ,  U/  is  the  acceleration  vector  field,  F“'  is  the 
external  force  applied  on  the  node  I  by  evaluating  the  corresponding  components  of  F  in 
Equation  (28),  and  F""  is  the  internal  force  vector  of  the  node  I .  Because  the  mass  matrix  is 
diagonal,  Equation  (30)  can  be  solved  without  factorizing  a  global  stiffness  matrix. 


To  bridge  the  FE  subregion  and  the  peridynamic  subregion,  we  introduce  an  interface  element.  A 
three-dimensional  interface  element  consisting  of  eight  FE  nodes  is  illustrated  in  Figure  5.  In  an 
interface  element,  a  number  of  peridynamic  nodes  are  embedded  for  the  calculation  of  coupling 
forces.  The  interacting  forces  between  embedded  peridynamic  nodes  and  peridynamic  nodes  out 
of  interface  elements  are  called  coupling  forces.  It  is  worth  noting  that  interactions  between 
embedded  peridynamic  nodes  are  not  considered  as  coupling  forces.  The  number  of  embedded 
peridynamic  nodes  is  determined  by  the  size  of  the  horizon,  and  there  should  be  sufficient 
embedded  nodes  within  the  horizon  of  nodes  near  the  interface  of  the  peridynamic  subregion  and 
the  FE  subregion  as  shown  in  Figure  4. 
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Figure  4:  Interface  element  for  the  coupling  of  FE  subregion  and  peridynamic  subregion. 


To  evaluate  coupling  forces,  each  embedded  peridynamic  node  represents  a  material  volume  of 
(Ax)3  inside  an  interface  element.  However,  embedded  peridynamic  nodes  are  not  involved  in 
the  global  equation.  In  other  words,  the  displacements  of  embedded  peridynamic  nodes  are  not 
calculated  by  solving  the  equation  of  motion.  Therefore,  the  mass  of  an  interface  element  is 
equally  distributed  to  FE  nodes  of  the  interface  element.  Consider  an  embedded  peridynamic 
node  subjected  to  the  coupling  force  fcp  as  shown  in  Figure  5,  the  coupling  force  is  then 
divided  to  FE  nodes  of  the  interface  element  by  means  of  shape  functions  as 

^'=^.(c,//,v/)fr,\  i  =  (31) 

where  $  is  the  shape  function  of  the  node  i  belonging  to  the  interface  element,  (i;,r/,y/)  are 

the  natural  coordinates  of  the  embedded  node  in  the  interface  element,  which  should  be 
determined  by  the  inverse  isoparametric  mapping.  We  designate  this  coupling  scheme  as  the 
VL-coupling  scheme  since  the  whole  volume  of  the  interface  element  is  subjected  to  coupling 
forces.  On  the  other  hand,  different  from  the  VL-coupling  scheme,  we  might  divide  coupling 
forces  only  to  the  FE  nodes  on  the  interface  segment  as  shown  in  Figure  6.  Therefore,  FE  nodes 
not  on  the  interface  segment  are  subjected  to  internal  forces  arising  from  the  element  stiffness 
only.  Since  the  interface  between  the  peridynamic  and  FE  subregions  is  similar  to  a  contact 
surface,  the  scheme  demonstrated  in  Figure  6  is  designated  as  the  CT-coupling  scheme.  To 
implement  the  CT-coupling  scheme,  interfaces  between  the  peridynamic  subregion  and  the  FE 
subregion  have  to  be  defined  prior  to  analyses.  Coupling  forces  on  embedded  nodes  are  divided 
to  those  FE  nodes  on  the  interface  segment  as  shown  in  Figure  6  by 

i7=^Vc)icp,  i  =  3, 4, 7, 8,  (32) 

where  $  is  the  shape  function  on  the  interface  segment,  and  (£ c,r/c )  are  the  natural 
coordinates  of  the  projection  of  an  embedded  node  onto  the  interface  segment. 

In  general,  the  equation  of  motion  for  FE  nodes  of  an  interface  element  is  written  as 

M,t,=  F/-tt  -F;"\  (33) 

where  F/xr  is  the  external  force  by  evaluating  Equation  (28),  and  the  internal  force  is  given  as 


(34) 
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in  which  [-]7  denotes  the  corresponding  components  of  a  vector  associated  with  the  node  I  of 

the  interface  element,  and  f)7'  is  the  summation  of  coupling  forces  on  the  node  /  .  The  explicit 

algorithm  is  employed  for  the  transient  dynamic  analyses.  Nodal  accelerations  are  calculated  first, 
and  nodal  velocities  and  displacements  are  updated  subsequently.  After  the  displacements  of  FE 
nodes  of  the  interface  elements  are  calculated,  the  displacements  of  embedded  peridynamic 
nodes  are  then  determined  by 

</,*')  U„  1  =  1,-,  8,  (35) 

in  which  are  the  natural  coordinates  of  an  embedded  peridynamic  node  in  the  interface 

element,  and  U(.  is  the  nodal  displacement  of  an  interface  element.  For  peridynamic  nodes  out 
of  the  interface  elements,  Equation  (10)  is  used  to  update  nodal  accelerations. 


4.2  Inverse  isoparametric  mapping 

To  couple  peridynamic  subregions  with  finite  element  subregions,  a  certain  number  of 
peridynamic  nodes  are  embedded  in  the  interface  elements.  If  the  Cartesian  coordinates  of  an 
embedded  peridynamic  node  are  known,  the  natural  coordinates  of  the  embedded  peridynamic 
node  in  an  interface  element  should  be  determined  by  the  inverse  isoparametric  mapping,  which 
is  essential  especially  for  random  discretizations. 


:  J  Interface  element 
^  FE  node 
0  PD  node 
0  embedded  PD  node 


Figure  5:  VL-coupling  scheme  that  divides  a  coupling  force  fcp  to  FE  nodes  of  the  interface  element. 


However,  the  inverse  isoparametric  mapping  from  the  Cartesian  coordinates  to  the  natural 
coordinates  is  nontrivial  since  equations  to  be  solved  are  nonlinear.  Murti  and  Valliappan  [62] 
presented  a  numerical  technique  by  bisecting  a  line  passing  a  point  and  a  node  of  known  natural 
coordinates,  and  this  method  was  extended  to  the  three-dimensional  space  by  Murti  et  al.  [63]. 
However,  the  bisection  method  has  its  limitations  [64], 


Figure  6:  CT-coupling  scheme  that  divides  a  coupling  force  fcp 


to  FE  nodes  on  the  interface  segment. 


A  more  generalized  approach  for  the  inverse  isoparametric  mapping  is  presented  by 
Chinnaswamy  et  al.  [64],  For  the  inverse  mapping  of  an  embedded  peridynamic  node  with 
known  Cartesian  coordinates  (x,  y,  z) ,  the  equation  to  be  solved  is  written  as 
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where  are  the  natural  coordinates  of  an  embedded  peridynamic  node  in  an  interface 

element  to  be  determined.  By  expanding  the  vector  f  in  Taylor’s  series  and  omitting  the  second 
and  higher  order  terms,  it  can  be  shown  that  [64] 
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where  I0  =  [<^0,70,^0]r  is  an  approximate  solution,  and  f0  is  the  vector  f  evaluated  at  the 
approximate  solution  I0 .  Equation  (37)  can  be  rewritten  as 
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and  it  can  be  simplified  as 
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The  updated  solution  of  I  is  used  as  the  value  of  I0  in  Equation  (39)  for  the  next  iteration.  A 
few  iterations  are  performed  till  the  solution  of  I  converges. 


If  the  CT-coupling  scheme  is  employed,  coupling  forces  on  embedded  peridynamic  nodes  are 
only  distributed  to  FE  nodes  on  the  interface  segment  as  shown  in  Figure  6.  Hence,  the  natural 
coordinates  of  the  projection  of  an  embedded  peridynamic  node  onto  the  interface  segment  have 
to  be  determined. 


Figure  7:  Projection  of  an  embedded  peridynamic  node  on  an  interface  segment. 

Let  t  be  the  position  vector  of  an  embedded  node  ns,  and  the  projection  of  ns  onto  the 
interface  segment  is  denoted  by  ns,  as  shown  in  Figure  7.  The  position  vector  r  of  the  point 
ns,  on  the  interface  segment  can  be  expressed  as 

r  =  +/3(4^c)e  3>  (41) 

where  (g c,r\c )  are  the  natural  coordinates  of  the  point  ns,  on  the  interface  segment,  and 

m,ri)  =  Yj>jxU  (42) 

y=i 

in  which  <f>.  is  the  shape  function  of  the  node  j  on  the  interface  segment.  The  natural 
coordinates  (£ c,rjc )  of  the  point  ns,  on  the  interface  segment  must  satisfy  [65] 

||-(t-r)  =  0,  (43) 

|L'(t-r)  =  0.  (44) 

orj 

However,  there  is  no  analytical  solution  to  Equations  (43)  and  (44).  To  solve  numerically,  a  few 
iterations  of  the  least-squares  projection  are  used  to  generate  an  initial  guess  as 
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A^A?7  =  -r, -(r-t)r  -(r-t), 


5+I=£+A£  17^=17,+ A 17.  (47) 

With  an  initial  guess,  Newton-Raphson  method  is  then  utilized  to  find  the  solution  of  Equations 
(43)  and  (44)  as  [66] 

HA<f;A  rj  =  -r ^  •  (r  -  t)r  7  •  (r  - 1),  (48) 
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4i+i  =  4i + 7i+i  =  ^ + A  n-  (50) 

The  solutions  of  <%M  and  rjM  are  used  to  update  the  value  of  the  position  vector  r ,  and  then 
Equation  (48)  is  evaluated  again  for  the  next  iteration.  The  converged  solutions  of  and  rjM 
are  the  natural  coordinates  of  the  point  ns,  on  the  interface  segment. 


5  Numerical  applications 

5. 1  One-dimensional  bar 

For  benchmarking,  the  present  coupling  approach  of  peridynamics  with  FEM  is  employed  to 
study  the  axial  deformation  of  a  one-dimensional  bar.  The  length  of  the  bar  is  9.5  mm,  and 
dimensions  of  the  cross  section  are  0.5  mm  by  0.5  mm.  Young’s  modulus  E  of  the  bar  is 
70  GPa,  and  the  density  is  2700  kg/m3.  Figure  8  shows  the  multiscale  discretization  of  the 
bar.  The  finite  element  mesh  size  is  1.5  mm,  and  the  conventional  bar  elements  are  utilized. 

EA 

The  stiffness  of  a  bar  is  k  =  -j- ,  where  A  is  the  cross-sectional  area  and  L  is  the  length  of 

the  bar.  The  peridynamic  grid  spacing  is  0.5  mm,  which  is  equal  to  the  width  of  the  bar.  By 
setting  the  horizon  to  2 Ax,  the  one-dimensional  micromodulus  q  is  5.6 xlO23  N/m6.  Two 
interface  elements  are  used  to  couple  peridynamic  and  FE  subregions,  and  each  interface  element 


has  two  embedded  peridynamic  nodes  for  the  calculation  of  coupling  forces  as  shown  in  Figure 
8. 
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Figure  8:  Discretization  of  a  one-dimensional  bar  for  coupling. 
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Figure  9:  Axial  displacement  along  the  bar  using  the  (a)  VL-coupling  scheme  and  (b)  CT-coupling 

scheme. 


A  tensile  loading  with  the  magnitude  of  F  =  175  N  is  applied  on  both  ends  of  the  bar  under  the 
quasi-static  condition.  The  force  is  gradually  increased  during  50,000  steps  with  the 

calculation  time  step  dt  =  5  x  10“8  s,  which  is  less  than  the  critical  time  step  for  the  explicit  time 
integration.  Figure  9(a)  shows  the  displacement  along  the  bar  using  the  VL-coupling  scheme. 
Nodal  displacements  in  the  peridynamic  subregion  show  good  agreement  with  the  quasi-static 
solution.  The  displacements  of  the  FE  nodes,  however,  show  small  discrepancies.  The  reason  is 
that  coupling  forces  in  an  interface  element  are  divided  to  all  FE  nodes  of  the  element  using  the 
VL-coupling  scheme  so  that  FE  nodes  at  the  interfaces  of  subregions  (x  =  ±1.75  mm)  only 
receive  partial  coupling  forces.  FE  nodes  at  the  other  end  of  the  interface  elements  (x  =  ±3.25 


mm)  receive  the  rest  of  coupling  forces,  and  are  also  subjected  to  internal  forces  arising  from  the 
element  stiffness.  Consequently,  the  displacements  of  FE  nodes  at  the  interfaces  are  slightly 
overestimated,  and  displacements  of  other  FE  nodes  are  underestimated  as  indicated  in  Figure 
9(a).  In  contrast,  the  solution  using  the  CT-coupling  scheme,  which  distributes  coupling  forces 
only  to  FE  nodes  at  interfaces,  is  almost  identical  to  the  quasi-static  solution  as  shown  in  Figure 
9(b).  Hence,  for  the  calculation  of  axial  displacement,  the  CT-coupling  scheme  is  more  effective 
than  the  VL-coupling  scheme  to  achieve  the  coupling  between  the  peridynamic  subregion  and 
the  FE  subregion. 


5.2  Three-dimensional  bar 


Interface  element  PD  subregion  Interface  element 

I*  1 1 1  ■  l  ■  l 


Gx 


Figure  10:  Three-dimensional  bar  subjected  to 
tension. 


A  three-dimensional  bar  subjected  to  tension  is 
examined  to  compare  the  solutions  of  the 
present  coupling  approach  and  the  classical 
(local)  elasticity  solutions.  The  dimensions  of 
the  bar  are  taken  to  be  10  mm  in  length,  7 
mm  in  width,  and  7  mm  in  thickness  as 
shown  in  Figure  10.  The  three-dimensional 
model  is  partitioned  into  two  FE  subregions 
and  one  peridynamic  subregion.  Each  FE 
subregion  consists  of  four  eight-node  solid 


interface  elements,  and  the  mesh  size  of  the  interface  element  is  3.5  mm.  The  peridynamic 
subregion  is  uniformly  discretized  with  the  grid  spacing  Ax  =  0.5  mm,  and  the  size  of  the 
horizon  is  set  to  <7  =  1.0  mm.  In  the  interface  element,  two  additional  layers  of  peridynamic 
nodes  are  embedded  along  the  longitudinal  direction  to  ensure  sufficient  nodes  in  the  horizon  of 
peridynamic  nodes  near  the  interfaces.  Tractions  at  both  end  surfaces  are  gradually  applied  up  to 
crx  =  700  MPa  during  100,000  steps  as  quasi-static  loading,  and  the  calculation  time  step  is  set 

to  dt  =  5xl0“8  s.  Young’s  modulus  of  the  bar  is  70  GPa,  Poisson’s  ratio  is  0.25,  and  the 
density  is  2700kg/m3. 

We  first  examine  the  numerical  solutions  using  the  VL-coupling  scheme.  The  longitudinal 
displacements  ux  at  three  measuring  positions  are  compared  to  the  quasi-static  results  as  shown 


in  Figure  11(a).  The  displacement  ux  at  the  position  px,  which  is  the  at  the  end  surface,  is 
underestimated  compared  to  the  quasi-static  solution.  On  the  other  hand,  the  displacement  at  p2 , 
which  is  at  the  interface  of  the  peridynamic  subregion  and  the  FE  subregion,  and  the 
displacement  at  p3 ,  which  is  inside  the  peridynamic  subregion,  show  good  agreement  with  the 
quasi-static  solutions.  To  look  into  Poisson  effect,  the  transverse  displacements  measured  at 
different  positions  are  plotted  in  Figure  11(b).  The  transverse  displacement  qx  at  the  end 
surface  is  smaller  than  the  quasi-static  value  since  the  longitudinal  displacement  at  the  end 
surface  is  underestimated.  Nevertheless,  the  transverse  displacement  q2  at  the  interface 
demonstrates  close  agreement  with  the  quasi-static  solution,  which  indicates  that  Poisson  effect 
is  preserved  at  the  interface.  The  transverse  displacement  at  q2 ,  which  is  inside  the  peridynamic 
subregion,  also  matches  the  quasi-static  solution. 

For  the  comparison  of  two  types  of  coupling  schemes,  the  CT-coupling  scheme,  which  achieves 
coupling  by  considering  the  interfaces  of  subregions  similar  to  contact  surfaces,  is  then  employed 
to  study  this  quasi-static  problem.  As  shown  in  Figure  12(a),  the  longitudinal  displacement  px 
at  the  end  surface  is  close  to  the  quasi-static  solution  with  the  error  less  than  2  %.  The 
longitudinal  displacements  at  interface  and  inside  the  peridynamic  subregion,  denoted  by  p2 

and  p2  respectively  in  Figure  12(a),  are  almost  identical  to  the  quasi-static  solutions.  The 

transverse  displacements  measured  at  different  positions  are  plotted  in  Figure  12(b).  With  the 
improved  result  in  the  longitudinal  displacement  at  the  end  surface,  the  transverse  displacement 
qx  at  the  end  surface  turns  to  be  very  close  to  the  quasi-static  solution.  On  the  other  hand,  the 

transverse  displacement  q2  at  the  interface  is  overestimated.  This  phenomenon  occurs  due  to 

the  reason  that  decompositions  of  coupling  forces  in  the  transverse  direction  are  divided  only  to 
FE  nodes  on  the  interface.  By  comparing  Figures  11  and  12,  it  is  observed  that  the  CT-coupling 
scheme  is  effective  in  resolving  displacements  normal  to  the  interface  of  peridynamic  and  FE 
subregions.  On  the  other  hand,  the  VL-coupling  scheme,  which  divides  coupling  forces  to  all  FE 
nodes  of  the  interface  element,  is  capable  to  preserve  Poisson  effect  at  the  interface. 
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Figure  1 1 :  Displacements  (a)  ux  and  (b)  u  at  different  positions  using  the  VL-coupling  scheme. 


Coordinates  of  the  measuring  position  are  given  in  parenthesis. 
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Figure  12:  Displacements  (a)  ux  and  (b)  u  at  different  positions  using  the  CT-coupling  scheme. 
Coordinates  of  the  measuring  position  are  given  in  parenthesis. 
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Figure  13:  Axial  displacement  along  the  edge  of  the  bar 
using  the  CT-coupling  scheme. 


Figure  13  shows  the  longitudinal 
displacement  along  the  edge  of  the  bar 
solved  by  the  CT-coupling  scheme.  It  is 
noted  that  a  smooth  curve  can  be 
obtained  if  the  nodal  displacements  are 
connected,  which  is  different  from  the 
result  in  [52]  where  a  jump  in 
displacement  is  observed  at  the 
interface.  The  strain  sx  distributions  in 

the  bar  and  at  the  interface  are  plotted  in 
Figure  14.  Nodal  strains  in  the  FE 
subregions  are  obtained  by  evaluating 


Equation  (23),  and  nodal  strains  in  the  peridynamic  subregion  are  calculated  as  the  average  bond 
stretch  within  the  grid  width  Ax  .  Considering  the  quasi-static  condition,  the  analytical  value  of 
strain  ex  is  0.01 .  The  strain  of  the  coupling  model  is  in  the  range  of  0.0099  to  0.0105  as 
shown  in  Figure  14,  which  agrees  with  the  quasi-static  solution. 
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Figure  14:  Strain  sx  distributions  (a)  in  the  bar  and  (b)  at  the  interface  when  the  applied  traction 

Gx  =700  MPa.  The  CT-coupling  scheme  is  used  in  the  simulation. 


The  essence  of  coupling  forces  is 
interactions  between  nodes  in  the 
peridynamic  subregion  and  material 
volumes  of  an  adjacent  continuous 
body  (i.e.  interface  elements) 
represented  by  embedded 
peridynamic  nodes.  Therefore,  the 
ratio  of  the  grid  spacing  of 
peridynamic  nodes  to  the  mesh  size 
of  interface  elements  should  not 
affect  the  results.  To  illustrate  this 
point,  Figure  15  shows  the 
summation  of  coupling  forces  in  the 
longitudinal  direction  on  the 
interface  elements  at  the  left  end  of  the  bar  as  the  applied  traction  increases.  As  indicated  by  the 
comparison  in  Figure  15,  differences  in  the  results  using  the  grid  spacing  Ax  =  1.0  mm  and 
Ax  =  0.5  mm  are  insignificant. 

5.3  Mixed  mode  fracture  in  a  tension-shear  specimen 

A  benchmark  problem  of  mixed  mode  crack  propagation  in  a  concrete  specimen  has  been 
investigated  experimentally  by  Nooru-Mohamed  et  at.  [19].  The  double-edge-notched  specimen 
is  illustrated  in  Figure  16(a).  The  dimensions  of  the  specimen  are  taken  to  be  200  mm  in  both 
length  and  height,  50  mm  in  width,  and  two  notches  at  edges  are  25  mm  in  length,  5  mm  in 
height,  and  50  mm  in  width.  For  the  numerical  study,  the  specimen  is  partitioned  into  two  FE 
subregions  and  one  peridynamic  subregion  as  shown  in  Figure  16(b).  The  FE  subregion  is 
discretized  with  two  mesh  sizes,  10  mmx  10  mmx  16.25  mm  and  10  mmx  10  nun  x  13  mm. 
The  peridynamic  subregion  is  discretized  with  the  grid  spacing  Ax  =  5.0  mm,  which  is  1/10 
of  the  specimen  thickness,  and  the  horizon  is  set  to  8  =  1.0  mm.  The  notches  in  the 
peridynamic  subregion  are  introduced  by  deleting  nodes  along  two  notches  and  removing  all 
bonds  across  the  notches.  Since  the  ratio  of  the  horizon  to  the  grid  spacing  is  equal  to  two,  two 


Figure  15:  Summation  of  coupling  forces  in  the  longitudinal 
direction  on  the  interface  elements  at  the  left  end  of  the  bar. 
The  CT-coupling  scheme  is  used  in  the  simulation. 


layers  of  peridynamic  nodes  adjacent  to  interfaces  are  embedded  in  interface  elements  for  the 
calculation  of  coupling  forces. 
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Figure  16:  Mixed  mode  fracture  test:  (a)  geometry  of  the  specimen;  (b)  subregions  of  the  coupling  model. 


Young’s  modulus,  Poisson’s  ratio,  and  fracture  energy  were  not  measured  in  the  experiments. 
Therefore,  we  adopt  the  material  properties  E  =  30  GPa  and  Gf  =110  J/m2  as  in  [10].  For 

Poisson’s  ratio,  it  is  estimated  as  v  =  0.2  in  the  numerical  studies  in  [10,  12],  and  v  =  0.3  is 
used  in  [67].  In  the  present  study,  Poisson’s  ratio  v  =  0.25  is  assumed.  The  density  is  calculated 
from  the  concrete  composition  given  in  [19]  as  p  =  2265  kg/m3 .  By  applying  Equation  (8),  the 
critical  bond  stretch  is  obtained  as  s0  =  5.5277  x  10  4 .  The  brittle  material  model  illustrated  in 

Figure  2  is  utilized  for  the  peridynamic  subregion,  and  the  linear  elastic  model  is  employed  for 
the  remaining  FE  subregions,  where  material  parameters  E  =  30  GPa  and  v  =  0.25  are  used. 


For  the  comparison,  the  load-path  lb  (specimen  46-05)  in  [19]  is  considered.  A  horizontal  shear 
force  of  10  kN  is  applied  first,  and  then  the  vertical  displacement  un  is  applied  on  the  top  and 
bottom  of  the  specimen  as  shown  in  Figure  16(a).  In  the  numerical  calculation,  the  time  step  is 
set  to  dt  =  lxlO-7  s,  which  is  less  than  the  critical  time  step  for  the  explicit  time  integration. 
The  vertical  displacement  un  is  applied  by  imposing  a  constant  velocity  of  10  mm/s  as 


quasi-static  loading.  We  first  examine  the  numerical  predictions  of  crack  paths  using  the 
VL-coupling  scheme.  The  cracks  initiate  near  the  notches  as  shown  in  Figure  17(a),  and 
propagate  along  the  horizontal  direction  for  about  50  mm.  As  the  boundary  displacement 
increases,  the  direction  of  propagation  changes  as  shown  in  Figures  17(b)  and  17(c).  At  the 
boundary  displacement  un  =  0.09  mm,  two  cracks  are  connected  as  shown  in  Figure  17(d).  The 

numerical  prediction  of  crack  paths  using  the  YL-coupling  scheme  shows  differences  with  the 
experimental  observations  in  [19]. 

We  then  apply  the  CT-coupling  scheme  for  the  numerical  simulation.  Interfaces  are  predefined  in 
the  reference  configuration,  and  there  are  280  interface  segments  with  the  given  mesh  sizes. 
Natural  coordinates  of  the  projections  of  embedded  nodes  onto  interface  segments  are  saved  to  a 
list.  In  the  subsequent  calculations,  the  natural  coordinates  of  projected  points  are  referred,  and 
the  coupling  forces  on  embedded  peridynamic  nodes  are  divided  to  FE  nodes  on  interface 
segments.  The  damage  evolution  is  shown  in  Figure  18.  Crack  initiation  occurs  at  the  left  and 
right  notches  as  shown  in  Figure  18(a),  and  materials  ahead  of  crack  tips  are  damaged  for  the 
length  of  around  15  mm  to  20  mm.  Due  to  the  angle  change  of  the  principal  stress,  cracks 
propagate  with  an  angle  as  shown  in  Figure  18(b).  As  the  boundary  displacement  increases,  two 
curvilinear  crack  paths  are  clearly  observed  in  Figure  18(c).  An  enclose  area  is  gradually  formed 
between  two  curvilinear  cracks  as  indicated  in  Figure  18(d).  The  numerical  prediction  of  crack 
paths  shows  agreement,  especially  for  the  lower  crack,  with  the  experimental  observation 
presented  in  [19],  which  is  illustrated  using  the  solid  line  in  Figure  18(d).  The  small 
discrepancies  appeared  in  Figure  18(d)  might  be  caused  by  the  perfectly  brittle  material  model 
adopted  in  the  numerical  simulation.  More  investigation  of  material  models  for  concrete  in  the 
framework  of  peridynamics  is  required  in  the  future.  Note  that  the  numerically  predicted  crack 
paths  are  in  perfect  symmetry.  On  the  contrary,  the  cracks  in  the  experiment  do  not  show  perfect 
symmetry,  and  we  might  speculate  on  symmetry  breaking  caused  by  the  slight  differences  in  the 
loading  condition  and  the  material  heterogeneity. 

Compared  to  the  numerical  results  using  the  YL-coupling  scheme,  the  results  using  the 
CT-coupling  scheme  demonstrate  better  agreement  with  the  experimental  observation.  The 


reason  for  the  inferiority  of  the  VL-coupling  scheme  in  this  mixed  mode  fracture  problem  is  that 
the  horizontal  displacement  at  edges  caused  by  the  constant  shear  loading  is  underestimated  after 
the  applied  vertical  displacement  at  boundaries  reaches  the  critical  value  for  the  crack  initiation 
at  two  notches.  Consequently,  the  opening-mode  fracture  dominates,  and  the  large  area  ahead  of 
crack  tips  is  damaged  in  the  plane  of  notches.  After  the  intact  region  in  the  middle  of  the 
specimen  becomes  relatively  small,  rotations  of  crack  paths  then  take  place  as  shown  in  Figure 
17. 
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Figure  17:  Numerical  prediction  of  crack  paths  using  the  VL-coupling  scheme.  Boundary  displacement: 
(a)  un  =  0.0220  mm;  (b)  un  —  0.0225  mm;  (c)  un  —  0.03  mm;  (d)  un  =  0.09  mm. 


(C) 


(d) 


Figure  18:  Numerical  prediction  of  crack  paths  using  the  CT-coupling  scheme.  Boundary  displacement:  (a) 
un  =  0.0220  mm;  (b)  un  =  0.0225  mm;  (c)  un  =  0.03  mm;  (d)  un  =  0.09  mm.  Solid  curves  are 
crack  paths  observed  in  the  experiment  [19]. 


6  Conclusions 

A  coupling  approach  of  discretized  peridynamics  with  FEM  is  developed  in  this  research.  To 
bridge  conventional  FE  subregions  and  peridynamic  subregions,  an  interface  element  is 
introduced.  The  proposed  coupling  approach  is  different  from  other  methods  in  the  sense  of 
direct  coupling  via  interface  elements.  Depending  on  the  size  of  the  horizon,  a  number  of 
peridynamic  nodes  are  embedded  in  an  interface  element.  The  embedded  peridynamic  nodes  are 
not  involved  in  the  global  equation,  but  essential  in  the  calculation  of  coupling  forces.  The 
coupling  forces  describe  interactions  between  embedded  peridynamic  nodes  in  interface 
elements  and  peridynamic  nodes  in  peridynamic  subregions.  Two  types  of  coupling  schemes  are 
introduced.  In  the  YL-coupling  scheme,  coupling  forces  on  embedded  peridynamic  nodes  are 
divided  to  FE  nodes  of  interface  elements.  On  the  other  hand,  coupling  forces  are  divided  only  to 
FE  nodes  on  interface  segments  in  the  CT-coupling  scheme.  The  inverse  isoparametric  mapping 
techniques  to  determined  the  natural  coordinates  of  embedded  peridynamic  nodes  in  the  interface 
elements  and  the  natural  coordinates  of  projected  points  on  the  interface  segments  are 
summarized. 

Numerical  simulations  are  conducted  to  compare  the  computational  results  using  the  coupling 
approach  to  the  classical  elasticity  solutions.  The  axial  deformation  of  a  one-dimensional  bar 
under  quasi-static  loading  is  studied.  It  is  found  that  the  displacements  at  interfaces  of  subregions 
are  slightly  overestimated,  and  the  displacements  in  the  FE  subregions  are  underestimated  using 
the  VL-coupling  scheme.  On  the  other  hand,  the  numerical  solution  using  the  CT-coupling 
scheme  is  almost  identical  to  the  quasi-static  solution  in  all  subregions  and  interfaces.  For 
three-dimensional  simulations,  a  bar  subjected  to  quasi-static  tension  is  partitioned  into  a 
peridynamic  subregion  and  two  FE  subregions  consisting  of  eight-node  interface  elements  at 
both  ends  of  the  bar.  By  measuring  displacements  at  different  positions,  the  CT-coupling  scheme 
is  found  to  be  effective  in  resolving  displacements  normal  to  the  interface  of  peridynamic  and  FE 
subregions,  and  the  YL-coupling  scheme  is  capable  to  preserve  Poisson  effect  at  the  interfaces. 
Longitudinal  strain  distributions  in  the  bar  and  at  the  interface  using  the  CT-coupling  scheme 
demonstrate  good  agreement  with  the  quasi-static  solution. 


The  last  numerical  example  is  the  mixed  mode  fracture  in  a  concrete  specimen  subjected  to 
quasi-static  loading.  The  region  where  failure  is  expected  is  modeled  using  peridynamics,  and 
the  remaining  region  is  modeled  using  conventional  FEM  to  reduce  the  computational  cost. 
Numerical  predictions  of  crack  paths  using  the  VL-coupling  scheme  and  the  CT-coupling 
scheme  are  studied.  Two  independent  curvilinear  crack  paths  are  observed  in  the  results  using 
the  CT-couping  scheme,  and  the  numerical  predictions  of  crack  patterns  are  close  to  the 
experimental  observations  presented  in  [19]. 
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